Importing Cell Profiler Images to Python

PI Sonja Aits, Lund University

Nikolay Oskolkov, SciLifeLab, NBIS Long Term Support, nikolay.oskolkov@scilifelab.se

Abstract

Here we will show how to convert Cell Profiler (Cellomics) format C01 to 8-bit tiff images and import those images into Python for further use with scikit-image and Convolutional Neural Networks (CNN). We will also demonstrate basic image processing steps with PIL / Pillow scikit-image Python modules.

Table of Contents:

Convert C01 to TIFF

Images which were received from Sonja Aits were in C01 format which is an intrinsic Cell Profiler (Cellomics) format.

Unfortunately, Cell Profiler is primarily a Windows program which is almost uninstallable under Linux. Following Cell Profiler instrictions here https://github.com/CellProfiler/CellProfiler/wiki/Source-Installation-(Ubuntu-14.04-LTS) to install it from source on Ubuntu 14.04 LTS which is my version, I ran into multiple problems with ilastik module https://www.ilastik.org/documentation/basics/installation.html which does not seem to be possible to install under Linux.

In [1]:
from IPython.display import Image
Image('/home/nikolay/WABI/S_Aits/Presentation/CellProfiler2Python/CellProfiler.png')
Out[1]:

Therefore I decided to use Bio-Formats library for processing biological microscopy images. More specifically, I used the command line tool bftools which can be downloaded here https://www.openmicroscopy.org/bio-formats/downloads/, bftools includes a very handy tool bfconvert which converts between over 100 formats supported by Bio-Formats, including C01 Cell Profiler format. To run bfconvert let us first list the C01 files and build names for corresponding tiff-files:

In [10]:
import glob
c01_files = glob.glob('/home/nikolay/WABI/S_Aits/co1_images/*.C01')
c01_files
Out[10]:
['/home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f01d1.C01',
 '/home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f00d1.C01',
 '/home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f01d0.C01',
 '/home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f00d0.C01',
 '/home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f01d2.C01',
 '/home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f00d2.C01']
In [18]:
tiff_files = [str(i.split('.C01')[0]) + '.tiff' for i in c01_files]
tiff_files = [i.replace('co1_images','tiff_images') for i in tiff_files]
tiff_files
Out[18]:
['/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d1.tiff',
 '/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d1.tiff',
 '/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d0.tiff',
 '/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff',
 '/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d2.tiff',
 '/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d2.tiff']

To run Unix commands such as bfconvert in Python environment we will use subprocess module:

In [27]:
import subprocess
for i,j in zip(c01_files, tiff_files):
    print('Converting {} to tiff format'.format(i))
    subprocess.run(['/home/nikolay/WABI/S_Aits/bftools/bfconvert', '-overwrite', '-nogroup', i, j])
Converting /home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f01d1.C01 to tiff format
Converting /home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f00d1.C01 to tiff format
Converting /home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f01d0.C01 to tiff format
Converting /home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f00d0.C01 to tiff format
Converting /home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f01d2.C01 to tiff format
Converting /home/nikolay/WABI/S_Aits/co1_images/MFGTMPcx7_170802000001_A01f00d2.C01 to tiff format

Now we have converted C01 to tiff-files, however the problem with the tiff-files is that they are 16-bit files which can not be open properly with many graph editors and viewers. We need to use the standard Linux convert command to generate 8-bit tiff-files:

In [57]:
%%bash
for i in $(ls /home/nikolay/WABI/S_Aits/tiff_images/*.tiff)
do
echo Conevrt file $i to 8-bit tiff format
convert $i -auto-level -depth 8 -define quantum:format=unsigned -type grayscale ${i::-5}_8bit.tiff
rm $i
mv ${i::-5}_8bit.tiff $i
done
Conevrt file /home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff to 8-bit tiff format
Conevrt file /home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d1.tiff to 8-bit tiff format
Conevrt file /home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d2.tiff to 8-bit tiff format
Conevrt file /home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d0.tiff to 8-bit tiff format
Conevrt file /home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d1.tiff to 8-bit tiff format
Conevrt file /home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d2.tiff to 8-bit tiff format

Hurray! Mission accomplished, we have converted C01-files to tiff format which can be opened by standard Windows / Mac / Linux graph viewers and editors. Now it is time to import those tiff-images into Python using special libraries for image editing.

Reading Images with ImageJ

Cell Profiler is hard to install on Linux but it is possible to install another Bio-Formats module which is called ImageJ. This way one can at least display the C01 images if not convert them to any other useful format. To display the images one needs to perform the following simple operations:

  • cd /home/nikolay/Downloads/ImageJ
  • ./ImageJ
  • Plugins -> Bio-Formats -> Bio-Formats Importer
  • Plugins -> Filters -> Color Transformer
In [1]:
from IPython.display import Image
Image('/home/nikolay/WABI/S_Aits/Presentation/CellProfiler2Python/ImageJ.png', width = 2000)
Out[1]:

Import Images to Python

Let us now try to import the tiff-images into Python for further image processing with Pillow and scikit-images Python modules. There are two common ways to import images into Python. First, we can use use scikit-image which immediately converts the image into numpy array:

In [81]:
from skimage import io
img = io.imread('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')[0][0][0]
In [82]:
img.shape
Out[82]:
(1104, 1104)
In [83]:
img
Out[83]:
array([[ 7,  9,  6, ...,  4,  4,  6],
       [ 9,  8,  7, ...,  6,  6,  5],
       [ 8, 10,  7, ...,  4,  6,  8],
       ..., 
       [32, 31, 29, ..., 35, 34, 39],
       [29, 31, 26, ..., 33, 36, 43],
       [22, 30, 26, ..., 30, 37, 39]], dtype=uint8)
In [98]:
import matplotlib.pyplot as plt
from skimage import io

fig = plt.figure(figsize=(20, 15))
for i,j in zip(tiff_files,range(len(tiff_files))):
    img = io.imread(i)[0][0][0]
    plt.subplot(231 + j)
    plt.imshow(img[:,:], cmap=plt.cm.gray)
    plt.title(i.split('/')[6].split('.')[0], fontsize = 20)
plt.tight_layout()
plt.show()

Superimposing Color Channels

Second, we can use PIL / Pillow module which is also very handy, but here we need to additionally convert the image into numpy array for future usage with Convolutional Neural Networks (CNNs). A great thing about PIL / Pillow module is that one can easily superimpose different color channels onto each other and get a colored image out of 3 black-and-white images. Let us demonstrate how to do it:

In [115]:
import matplotlib.pyplot as plt
from PIL import Image
im1 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
im2 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d1.tiff')
im3 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d2.tiff')

fig = plt.figure(figsize=(20, 15))
plt.figure(figsize = (12, 8))
plt.subplot(131)
plt.imshow(im1)
plt.subplot(132)
plt.imshow(im2)
plt.subplot(133)
plt.imshow(im3)
plt.tight_layout()
plt.show()
<Figure size 1440x1080 with 0 Axes>

Again, we can easily convert images into numpy arrays:

In [110]:
import numpy as np
imarray1 = np.array(im1)
imarray2 = np.array(im2)
imarray3 = np.array(im3)
imarray1.shape
Out[110]:
(1104, 1104)
In [111]:
imarray1
Out[111]:
array([[ 7,  9,  6, ...,  4,  4,  6],
       [ 9,  8,  7, ...,  6,  6,  5],
       [ 8, 10,  7, ...,  4,  6,  8],
       ..., 
       [32, 31, 29, ..., 35, 34, 39],
       [29, 31, 26, ..., 33, 36, 43],
       [22, 30, 26, ..., 30, 37, 39]], dtype=uint8)

And finally, let us merge black-and-white images from R, G and B color channels and display the result which will be a colored Image.

In [112]:
newImage = Image.merge("RGB", [im1,im2,im3])
In [113]:
np.array(newImage).shape
Out[113]:
(1104, 1104, 3)
In [114]:
plt.figure(figsize = (12, 8))
plt.imshow(newImage)
plt.show()

Image Processing

Image processing is the process of creating a new image from an existing image, typically simplifying or enhancing the content in some way. It is a type of digital signal processing and is not concerned with understanding the content of an image. Examples of image processing include:

  • Normalizing photometric properties of the image, such as brightness or color.
  • Cropping the bounds of the image, such as centering an object in a photograph.
  • Removing digital noise from an image, such as digital artifacts from low light levels.

We will perform image processing using Pillow Python module. Pillow is a PIL library that supports Python 3 and is the preferred modern library for image manipulation in Python. It is even required for simple image loading and saving in other Python scientific libraries such as SciPy and Matplotlib.

In [48]:
from PIL import Image
test_image = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
print(test_image.format)
print(test_image.mode)
print(test_image.size)
TIFF
L
(1104, 1104)
In [49]:
import matplotlib.pyplot as plt
plt.imshow(test_image)
plt.show()

Here mode 'L' means that this is an RGB image which was converted to grayscale. Another way to do image processing is in Matplotlib as well but it will run Pillow unders the hood.

In [50]:
# load and display an image with Matplotlib
from matplotlib import image
from matplotlib import pyplot
# load image as pixel array
data = image.imread('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
# summarize shape of the pixel array
print(data.dtype)
print(data.shape)
# display the array of pixels as an image
pyplot.imshow(data)
pyplot.show()
uint8
(1104, 1104)

We can access the pixel data from a Pillow Image. Perhaps the simplest way is to construct a NumPy array and pass in the Image object. The process can be reversed converting a given array of pixel data into a Pillow Image object using the Image.fromarray() function. This can be useful if image data is manipulated as a NumPy array and you then want to save it later as a PNG or JPEG file.

In [36]:
from numpy import asarray
data = asarray(test_image)
# summarize shape
print(data.shape)
# create Pillow image
image2 = Image.fromarray(data)
# summarize image details
print(image2.format)
print(image2.mode)
print(image2.size)
(1104, 1104)
None
L
(1104, 1104)

Both approaches are effective for loading image data into NumPy arrays, although the Matplotlib imread() function uses fewer lines of code than loading and converting a Pillow Image object and may be preferred. For example, you could easily load all images in a directory as a list as follows:

In [38]:
# load all images in a directory
from os import listdir
from matplotlib import image
import matplotlib.pyplot as plt
# load all images in a directory
loaded_images = list()
fig = plt.figure(figsize=(20, 15))
tiff_files = listdir('/home/nikolay/WABI/S_Aits/tiff_images/')
for filename,j in zip(tiff_files,range(len(tiff_files))):
    # load image
    img_data = image.imread('/home/nikolay/WABI/S_Aits/tiff_images/' + filename)
    plt.subplot(231 + j)
    plt.imshow(img_data)
    # store loaded image
    loaded_images.append(img_data)
    print('> loaded %s %s' % (filename, img_data.shape))
plt.show()
> loaded MFGTMPcx7_170802000001_A01f00d0.tiff (1104, 1104)
> loaded MFGTMPcx7_170802000001_A01f01d1.tiff (1104, 1104)
> loaded MFGTMPcx7_170802000001_A01f00d1.tiff (1104, 1104)
> loaded MFGTMPcx7_170802000001_A01f01d2.tiff (1104, 1104)
> loaded MFGTMPcx7_170802000001_A01f00d2.tiff (1104, 1104)
> loaded MFGTMPcx7_170802000001_A01f01d0.tiff (1104, 1104)

An image object can be saved by calling the save() function. This can be useful if you want to save an image in a different format, in which case the ‘format‘ argument can be specified, such as PNG, GIF, or PEG. Saving images is useful if you perform some data preparation on the image before modeling. One example is converting color images (RGB channels) to grayscale (1 channel). There are a number of ways to convert an image to grayscale, but Pillow provides the convert() function and the mode ‘L‘ will convert an image to grayscale.

In [43]:
print(test_image.mode)
gs_image = test_image.convert(mode='L')
print(gs_image.mode)
np.array(gs_image).shape
L
L
Out[43]:
(1104, 1104)
In [46]:
plt.imshow(test_image)
plt.show()
In [45]:
plt.imshow(gs_image)
plt.show()

In our case, the images were already in the grayscale mode, so conversion did not change anything, i.e. mode was 'L' before conversion and it remained 'L' after conversion.

Our images are of a very high resolution which might be a problem for the CNN as those images demand a lot of RAM even if one uses minibatches. Here we will show how to downsize all the images and save them with a lower resolution. For example, the test photograph we have been working with has the width and height of (1104, 1104). We can resize it to (100, 100), in which case the largest dimension, e.g. the width, will be reduced to 100, and the height will be scaled in order to retain the aspect ratio of the image. In our case, we have equal width and height so the final dimension will be (100, 100).

In [53]:
import os
import numpy as np
from PIL import Image
tiff_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_images/')
tiff_downsized_files = os.mkdir('/home/nikolay/WABI/S_Aits/tiff_downsized_images/')
for filename in tiff_files:
    img_data = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/' + filename)
    print('> loaded %s %s' % (filename, np.array(img_data).shape))
    img_data.thumbnail((100,100))
    print('> downsized image: {}'.format(np.array(img_data).shape))
    print('> saving downsized image {}'.format(filename))
    img_data.save('/home/nikolay/WABI/S_Aits/tiff_downsized_images/' + filename)
> loaded MFGTMPcx7_170802000001_A01f00d0.tiff (1104, 1104)
> downsized image: (100, 100)
saving downsized image MFGTMPcx7_170802000001_A01f00d0.tiff
> loaded MFGTMPcx7_170802000001_A01f01d1.tiff (1104, 1104)
> downsized image: (100, 100)
saving downsized image MFGTMPcx7_170802000001_A01f01d1.tiff
> loaded MFGTMPcx7_170802000001_A01f00d1.tiff (1104, 1104)
> downsized image: (100, 100)
saving downsized image MFGTMPcx7_170802000001_A01f00d1.tiff
> loaded MFGTMPcx7_170802000001_A01f01d2.tiff (1104, 1104)
> downsized image: (100, 100)
saving downsized image MFGTMPcx7_170802000001_A01f01d2.tiff
> loaded MFGTMPcx7_170802000001_A01f00d2.tiff (1104, 1104)
> downsized image: (100, 100)
saving downsized image MFGTMPcx7_170802000001_A01f00d2.tiff
> loaded MFGTMPcx7_170802000001_A01f01d0.tiff (1104, 1104)
> downsized image: (100, 100)
saving downsized image MFGTMPcx7_170802000001_A01f01d0.tiff

Simple image manipulation can be used to create new versions of images that, in turn, can provide a richer training dataset when modeling. Generally, this is referred to as data augmentation and may involve creating flipped, rotated, cropped, or other modified versions of the original images with the hope that the algorithm will learn to extract the same features from the image data regardless of where they might appear.

In [57]:
# horizontal flip
hoz_flip = test_image.transpose(Image.FLIP_LEFT_RIGHT)
# vertical flip
ver_flip = test_image.transpose(Image.FLIP_TOP_BOTTOM)
# plot all three images using matplotlib
fig = plt.figure(figsize=(20, 15))
pyplot.subplot(131)
pyplot.title('Original Image', fontsize=20)
pyplot.imshow(test_image)
pyplot.subplot(132)
pyplot.imshow(hoz_flip)
pyplot.title('Horizontal Flip', fontsize=20)
pyplot.subplot(133)
pyplot.imshow(ver_flip)
pyplot.title('Vertical Flip', fontsize=20)
pyplot.show()

Let us now create a simple augmented data set byt flipping each of the 6 downsized images horizontally and vertically:

In [58]:
import os
import numpy as np
from PIL import Image
tiff_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_downsized_images/')
tiff_augmented_files = os.mkdir('/home/nikolay/WABI/S_Aits/tiff_augmented_images/')
for filename in tiff_files:
    img_data = Image.open('/home/nikolay/WABI/S_Aits/tiff_downsized_images/' + filename)
    print('> loaded %s %s' % (filename, np.array(img_data).shape))
    hoz_flip = img_data.transpose(Image.FLIP_LEFT_RIGHT)
    ver_flip = img_data.transpose(Image.FLIP_TOP_BOTTOM)
    print('> saving augmented images for {}'.format(filename))
    img_data.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/' + filename)
    hoz_flip.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/aug_hoz_' + filename)
    ver_flip.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/aug_ver_' + filename)
> loaded MFGTMPcx7_170802000001_A01f00d0.tiff (100, 100)
> saving augmented images for MFGTMPcx7_170802000001_A01f00d0.tiff
> loaded MFGTMPcx7_170802000001_A01f01d1.tiff (100, 100)
> saving augmented images for MFGTMPcx7_170802000001_A01f01d1.tiff
> loaded MFGTMPcx7_170802000001_A01f00d1.tiff (100, 100)
> saving augmented images for MFGTMPcx7_170802000001_A01f00d1.tiff
> loaded MFGTMPcx7_170802000001_A01f01d2.tiff (100, 100)
> saving augmented images for MFGTMPcx7_170802000001_A01f01d2.tiff
> loaded MFGTMPcx7_170802000001_A01f00d2.tiff (100, 100)
> saving augmented images for MFGTMPcx7_170802000001_A01f00d2.tiff
> loaded MFGTMPcx7_170802000001_A01f01d0.tiff (100, 100)
> saving augmented images for MFGTMPcx7_170802000001_A01f01d0.tiff

An image can be rotated using the rotate() function and passing in the angle for the rotation. The function offers additional control such as whether or not to expand the dimensions of the image to fit the rotated pixel values (default is to clip to the same size), where to center the rotation the image (default is the center), and the fill color for pixels outside of the image (default is black).

In [59]:
# create rotated versions of an image
from PIL import Image
from matplotlib import pyplot
fig = plt.figure(figsize=(20, 15))
pyplot.subplot(131)
pyplot.title('Original Image', fontsize=20)
pyplot.imshow(test_image)
pyplot.subplot(132)
pyplot.imshow(test_image.rotate(45))
pyplot.title('Rotate 45 Degrees', fontsize=20)
pyplot.subplot(133)
pyplot.imshow(test_image.rotate(90))
pyplot.title('Rotate 90 Degrees', fontsize=20)
pyplot.show()

Let us now make a funny thing. We will create an augmented data set by rotating each image from 0 to 360 degrees by 2 degree step. Therefore each image will be replicated by 180 augmented images. Since we have essentially two photographs, 3 channels each, the first photograph will be represented by 180x3 = 540 images, and the second photograph will be represented by other 540 images. The total data set thus will be 1080 images, where half comes from first photograph and the other half comes from the second photograph. Later we can mark augmented images originated from the first photograph as "cases" and the images from the second photograph as "controls" and run a simple case-control image classification to see if the CNN is capable of discriminating between the two photographs. But firss tlet us create the 1080 imges augmented data set:

In [70]:
import os
import numpy as np
from PIL import Image
tiff_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_downsized_images/')
tiff_augmented_files = os.mkdir('/home/nikolay/WABI/S_Aits/tiff_augmented_images/')
for filename in tiff_files:
    print('Augmenting image {}'.format(filename))
    for i in range(0,360,2):
        img_data = Image.open('/home/nikolay/WABI/S_Aits/tiff_downsized_images/' + filename)
        rot_image = img_data.rotate(i)
        rot_image.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/aug_' + str(i) + '_' + filename)
Augmenting image MFGTMPcx7_170802000001_A01f00d0.tiff
Augmenting image MFGTMPcx7_170802000001_A01f01d1.tiff
Augmenting image MFGTMPcx7_170802000001_A01f00d1.tiff
Augmenting image MFGTMPcx7_170802000001_A01f01d2.tiff
Augmenting image MFGTMPcx7_170802000001_A01f00d2.tiff
Augmenting image MFGTMPcx7_170802000001_A01f01d0.tiff
In [71]:
%%bash
cd /home/nikolay/WABI/S_Aits/tiff_augmented_images/
ls *.tiff | wc -l
1080

We can see that we got 1080 augmented images.

In [72]:
%%bash
cd /home/nikolay/WABI/S_Aits/tiff_augmented_images/
find . -name '*f00*' | wc -l
find . -name '*f01*' | wc -l
540
540

Out of them 540 have the 'f00' suffix meaning they belong to the first photograph, and 540 have 'f01' suffix implying they belong to the second photograph. This will be our case-control data set which we will split into trainig and test data sets.

What is CNN and how it works?

Convolutional Neural Network (CNN) is a type of an Artificial Neural Network (ANN) which makes sense to use when dat have a spatial structure. This can be a matrix where columns and rows are ordered and their positions can not be changed, i.e. resampling (permutation, bootstrapping) does not make sense. An image is an example of such an ordered matrix since each pixel has a well-defined neighbour along x- and y-axes. In principle, a text, speach or any Time Series experiment can also be viewed as an ordered matrix but it is a less common practice to use CNN for those types of data. There are three types of layers in a Convolutional Neural Network:

  1. Convolutional Layers.
  2. Pooling Layers.
  3. Fully-Connected Layers.

Fully-Connected layer is a standard Dense layer like in Multilayer Perceptron (MLP), in contrast the other two types of layers are CNN-specific. Both Convolutional and Pooling layers have neurons which are called filters, these filters make some transformation of the input data. The filters can "see" only a limited part of the image, e.g. 5x5 pixel region which is called receptive field. One can think about a receptive field as an input for a neuron. The receptive field window can be "dragged" along horizontal or vertical axis of the image with a stride step and produce multiple receptive fileds. After each filter / neuron was applied to the receptive fields, the output of each neuron will be an ordered matrix which is called a feature map. This terminology is applicable for both Convolutional and Pooling layers, the difference between them is in how texactly the filters / neurons do transformation of the data.

In [1]:
from IPython.display import Image
Image('/home/nikolay/WABI/S_Aits/Presentation/ObjectDetection/CNN.png', width = 2000)
Out[1]:

Let us go through those layers one-by-one and explain how each of them works, on the way we will introduce common terminology used in CNNs.

Convolutional Layers

Suppose the input dimensions of the image are (1, 100, 100), where the picture has 100x100 pixels resolution and one color channel with values varying from 0 to 255, i.e. it is a grayscale image. Let us apply 100 filters to the image with a 3x3 receptive filed and stride 1, meaning we shift the receptive field window 1 horizontal or vertical position. Each filter applies a linear combintation to the values of the receptive field, like in the figure above. The weights of the filters are learned by minimizing the cost function via back-propogation like in normal ANNs. The output of each neuron of the Convolutional Layer is a feature map with "treansformed" color values for each pixel of the original image, so in total we will have 100 feature maps, each of them has slightly reduced dimensions (something like 95x95) compared to the original image.

Pooling Layers

Pooling Layers make sense to apply to the feature maps (output) of the Convolutional Layer. A receptive filed of a Pooling Layer typicaly has dimensions 2x2 applied with a stride of 2 and the filter / neuron simply takes the maximum value out of the 4 digits in the 2x2 receptive filed. The number of neurons / filters in the Pooling Layer is always equal to the number of neurons in the Convolutional Layer. As a result, Pooling Layer produces 100 feature maps, each of them has twice reduced dimension, i.e. 50x50, compared to the original image.

Fully-Connected Layers

Finally, after applying Covolutional and Pooling layers sequntially a couple of times, the output is flattened, i.e. reshaped into a one-dimensional array and input into an output layer in a traditional ANN fully-connected manner. The output layer generates predictions and has a softmax or sigmoid (for binary classification problems) activation functions applied to each neuron.

CIFAR10 Image Classification with CNN

Here for warming up purposes we will use a Convolutional Neural Network (CNN) built using Keras for image classification task using CIFAR10 image data set. The CIFAR-10 dataset consists of 60,000 photos divided into 10 classes (hence the name CIFAR-10). Classes include common objects such as airplanes, automobiles, birds, cats and so on. The dataset is split in a standard way, where 50,000 images are used for training a model and the remaining 10,000 for evaluating its performance. The photos are in color with red, green and blue components, but are small measuring 32 by 32 pixel squares.

In [118]:
# Plot ad hoc CIFAR10 instances
from keras.datasets import cifar10
from matplotlib import pyplot
from scipy.misc import toimage
# load data
(X_train, y_train), (X_test, y_test) = cifar10.load_data()
# create a grid of 3x3 images
for i in range(0, 9):
    pyplot.subplot(330 + 1 + i)
    pyplot.imshow(toimage(X_train[i]))
# show the plot
pyplot.show()
Using TensorFlow backend.
Downloading data from https://www.cs.toronto.edu/~kriz/cifar-10-python.tar.gz
170500096/170498071 [==============================] - 33s 0us/step
In [121]:
# Simple CNN model for CIFAR-10
import numpy
from keras.datasets import cifar10
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.layers import Flatten
from keras.constraints import maxnorm
from keras.optimizers import SGD
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.utils import np_utils
from keras import backend as K
K.set_image_dim_ordering('th')

# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)

# load data
(X_train, y_train), (X_test, y_test) = cifar10.load_data()

# normalize inputs from 0-255 to 0.0-1.0
X_train = X_train.astype('float32')
X_test = X_test.astype('float32')
X_train = X_train / 255.0
X_test = X_test / 255.0

# one hot encode outputs
y_train = np_utils.to_categorical(y_train)
y_test = np_utils.to_categorical(y_test)
num_classes = y_test.shape[1]

# Create the model
model = Sequential()
model.add(Conv2D(32, (3, 3), input_shape=(3, 32, 32), padding='same', activation='relu', 
                 kernel_constraint=maxnorm(3)))
model.add(Dropout(0.2))
model.add(Conv2D(32, (3, 3), activation='relu', padding='same', kernel_constraint=maxnorm(3)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(512, activation='relu', kernel_constraint=maxnorm(3)))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='softmax'))
# Compile model
epochs = 25
lrate = 0.01
decay = lrate/epochs
sgd = SGD(lr=lrate, momentum=0.9, decay=decay, nesterov=False)
model.compile(loss='categorical_crossentropy', optimizer=sgd, metrics=['accuracy'])
print(model.summary())
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_1 (Conv2D)            (None, 32, 32, 32)        896       
_________________________________________________________________
dropout_1 (Dropout)          (None, 32, 32, 32)        0         
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 32, 32, 32)        9248      
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 32, 16, 16)        0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 8192)              0         
_________________________________________________________________
dense_1 (Dense)              (None, 512)               4194816   
_________________________________________________________________
dropout_2 (Dropout)          (None, 512)               0         
_________________________________________________________________
dense_2 (Dense)              (None, 10)                5130      
=================================================================
Total params: 4,210,090
Trainable params: 4,210,090
Non-trainable params: 0
_________________________________________________________________
None
In [122]:
# Fit the model
model.fit(X_train, y_train, validation_data=(X_test, y_test), epochs=epochs, batch_size=32)
# Final evaluation of the model
scores = model.evaluate(X_test, y_test, verbose=0)
print("Accuracy: %.2f%%" % (scores[1]*100))
Train on 50000 samples, validate on 10000 samples
Epoch 1/25
50000/50000 [==============================] - 673s 13ms/step - loss: 1.7131 - acc: 0.3824 - val_loss: 1.3476 - val_acc: 0.5166
Epoch 2/25
50000/50000 [==============================] - 680s 14ms/step - loss: 1.3214 - acc: 0.5249 - val_loss: 1.1938 - val_acc: 0.5752
Epoch 3/25
50000/50000 [==============================] - 726s 15ms/step - loss: 1.1457 - acc: 0.5933 - val_loss: 1.0786 - val_acc: 0.6140
Epoch 4/25
50000/50000 [==============================] - 772s 15ms/step - loss: 1.0225 - acc: 0.6380 - val_loss: 1.0171 - val_acc: 0.6418
Epoch 5/25
50000/50000 [==============================] - 746s 15ms/step - loss: 0.9190 - acc: 0.6742 - val_loss: 0.9699 - val_acc: 0.6567
Epoch 6/25
50000/50000 [==============================] - 939s 19ms/step - loss: 0.8339 - acc: 0.7056 - val_loss: 0.9516 - val_acc: 0.6620
Epoch 7/25
50000/50000 [==============================] - 652s 13ms/step - loss: 0.7630 - acc: 0.7326 - val_loss: 0.9101 - val_acc: 0.6792
Epoch 8/25
50000/50000 [==============================] - 707s 14ms/step - loss: 0.6963 - acc: 0.7523 - val_loss: 0.9014 - val_acc: 0.6853
Epoch 9/25
50000/50000 [==============================] - 689s 14ms/step - loss: 0.6376 - acc: 0.7740 - val_loss: 0.9053 - val_acc: 0.6905
Epoch 10/25
50000/50000 [==============================] - 698s 14ms/step - loss: 0.5820 - acc: 0.7945 - val_loss: 0.8981 - val_acc: 0.6950
Epoch 11/25
50000/50000 [==============================] - 690s 14ms/step - loss: 0.5373 - acc: 0.8102 - val_loss: 0.9207 - val_acc: 0.6923
Epoch 12/25
50000/50000 [==============================] - 673s 13ms/step - loss: 0.5026 - acc: 0.8228 - val_loss: 0.9116 - val_acc: 0.6955
Epoch 13/25
 4192/50000 [=>............................] - ETA: 9:27 - loss: 0.4547 - acc: 0.8435
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-122-db6b14e9fa97> in <module>()
      1 # Fit the model
----> 2 model.fit(X_train, y_train, validation_data=(X_test, y_test), epochs=epochs, batch_size=32)
      3 # Final evaluation of the model
      4 scores = model.evaluate(X_test, y_test, verbose=0)
      5 print("Accuracy: %.2f%%" % (scores[1]*100))

~/miniconda3/lib/python3.6/site-packages/keras/models.py in fit(self, x, y, batch_size, epochs, verbose, callbacks, validation_split, validation_data, shuffle, class_weight, sample_weight, initial_epoch, **kwargs)
    891                               class_weight=class_weight,
    892                               sample_weight=sample_weight,
--> 893                               initial_epoch=initial_epoch)
    894 
    895     def evaluate(self, x, y, batch_size=32, verbose=1,

~/miniconda3/lib/python3.6/site-packages/keras/engine/training.py in fit(self, x, y, batch_size, epochs, verbose, callbacks, validation_split, validation_data, shuffle, class_weight, sample_weight, initial_epoch, steps_per_epoch, validation_steps, **kwargs)
   1629                               initial_epoch=initial_epoch,
   1630                               steps_per_epoch=steps_per_epoch,
-> 1631                               validation_steps=validation_steps)
   1632 
   1633     def evaluate(self, x=None, y=None,

~/miniconda3/lib/python3.6/site-packages/keras/engine/training.py in _fit_loop(self, f, ins, out_labels, batch_size, epochs, verbose, callbacks, val_f, val_ins, shuffle, callback_metrics, initial_epoch, steps_per_epoch, validation_steps)
   1211                     batch_logs['size'] = len(batch_ids)
   1212                     callbacks.on_batch_begin(batch_index, batch_logs)
-> 1213                     outs = f(ins_batch)
   1214                     if not isinstance(outs, list):
   1215                         outs = [outs]

~/miniconda3/lib/python3.6/site-packages/keras/backend/tensorflow_backend.py in __call__(self, inputs)
   2330         updated = session.run(self.outputs + [self.updates_op],
   2331                               feed_dict=feed_dict,
-> 2332                               **self.session_kwargs)
   2333         return updated[:len(self.outputs)]
   2334 

~/miniconda3/lib/python3.6/site-packages/tensorflow/python/client/session.py in run(self, fetches, feed_dict, options, run_metadata)
    893     try:
    894       result = self._run(None, fetches, feed_dict, options_ptr,
--> 895                          run_metadata_ptr)
    896       if run_metadata:
    897         proto_data = tf_session.TF_GetBuffer(run_metadata_ptr)

~/miniconda3/lib/python3.6/site-packages/tensorflow/python/client/session.py in _run(self, handle, fetches, feed_dict, options, run_metadata)
   1122     if final_fetches or final_targets or (handle and feed_dict_tensor):
   1123       results = self._do_run(handle, final_targets, final_fetches,
-> 1124                              feed_dict_tensor, options, run_metadata)
   1125     else:
   1126       results = []

~/miniconda3/lib/python3.6/site-packages/tensorflow/python/client/session.py in _do_run(self, handle, target_list, fetch_list, feed_dict, options, run_metadata)
   1319     if handle is None:
   1320       return self._do_call(_run_fn, self._session, feeds, fetches, targets,
-> 1321                            options, run_metadata)
   1322     else:
   1323       return self._do_call(_prun_fn, self._session, handle, feeds, fetches)

~/miniconda3/lib/python3.6/site-packages/tensorflow/python/client/session.py in _do_call(self, fn, *args)
   1325   def _do_call(self, fn, *args):
   1326     try:
-> 1327       return fn(*args)
   1328     except errors.OpError as e:
   1329       message = compat.as_text(e.message)

~/miniconda3/lib/python3.6/site-packages/tensorflow/python/client/session.py in _run_fn(session, feed_dict, fetch_list, target_list, options, run_metadata)
   1304           return tf_session.TF_Run(session, options,
   1305                                    feed_dict, fetch_list, target_list,
-> 1306                                    status, run_metadata)
   1307 
   1308     def _prun_fn(session, handle, feed_dict, fetch_list):

KeyboardInterrupt: 
In [138]:
X_train.shape
Out[138]:
(50000, 3, 32, 32)
In [137]:
y_train.shape
Out[137]:
(50000, 10)

Cell Image Classification with CNN

Finally here we will implement a vanilla CNN for classification of biological cell images. First, we need read the images using the Pillow Python module and convert them to numpy arrays. Important to note that we will consider all the 6 images as independent even though it is in fact 2 independent images that have 3 RGB channels each:

In [8]:
from PIL import Image
im1 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
im2 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d1.tiff')
im3 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d2.tiff')
im4 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d0.tiff')
im5 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d1.tiff')
im6 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d2.tiff')
#im1 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f00d0.png')
#im2 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f00d1.png')
#im3 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f00d2.png')
#im4 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f01d0.png')
#im5 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f01d1.png')
#im6 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f01d2.png')
In [9]:
import numpy as np
imarray1 = np.array(im1)
imarray2 = np.array(im2)
imarray3 = np.array(im3)
imarray4 = np.array(im4)
imarray5 = np.array(im5)
imarray6 = np.array(im6)

Looking at the dimensions of the arrays we conclude that each of them is a 4-dimensional array. First dimension corresponds to the number of samples, which is 1 for each image; second corresponds to the channel, which is 1 since we treat R, G and B channels as independent; finally, the last two dimensions correspond to the numbers of pixels in each 2D image.

In [10]:
imarray1_f = np.expand_dims(np.expand_dims(imarray1, axis=0), axis=0)
imarray2_f = np.expand_dims(np.expand_dims(imarray2, axis=0), axis=0)
imarray3_f = np.expand_dims(np.expand_dims(imarray3, axis=0), axis=0)
imarray4_f = np.expand_dims(np.expand_dims(imarray4, axis=0), axis=0)
imarray5_f = np.expand_dims(np.expand_dims(imarray5, axis=0), axis=0)
imarray6_f = np.expand_dims(np.expand_dims(imarray6, axis=0), axis=0)
imarray6_f.shape
Out[10]:
(1, 1, 1104, 1104)

Here we concatenate 4 images and build a training data set, while the rest two concatenated images will be for testing.

In [11]:
X_image_train = np.concatenate((imarray1_f, imarray2_f, imarray3_f, imarray4_f), axis = 0)
X_image_test = np.concatenate((imarray5_f, imarray6_f), axis = 0)
X_image_train.shape
Out[11]:
(4, 1, 1104, 1104)

We assign random binary classes, i.e. 0 and 1, to the training and test data sets.

In [12]:
y_image_train = np.array([0, 0, 1, 1])
y_image_test = np.array([0, 1])

Finally, we will construct the model and normalize the data:

In [13]:
import numpy
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.layers import Flatten
from keras.constraints import maxnorm
from keras.optimizers import SGD
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.utils import np_utils
from keras import backend as K
K.set_image_dim_ordering('th')

# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)

# normalize inputs from 0-255 to 0.0-1.0
X_image_train = X_image_train.astype('float32')
X_image_test = X_image_test.astype('float32')
X_image_train = X_image_train / 255.0
X_image_test = X_image_test / 255.0

# one hot encode outputs
y_image_train = np_utils.to_categorical(y_image_train)
y_image_test = np_utils.to_categorical(y_image_test)
num_classes = y_image_test.shape[1]
Using TensorFlow backend.
In [17]:
# Create the model
model = Sequential()
model.add(Conv2D(100, (3, 3), input_shape=(1, 1104, 1104), padding='same', activation='relu', 
                 kernel_constraint=maxnorm(3)))
model.add(Dropout(0.2))
model.add(Conv2D(100, (3, 3), activation='relu', padding='same', kernel_constraint=maxnorm(3)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(10, activation='relu', kernel_constraint=maxnorm(3)))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='sigmoid'))
# Compile model
epochs = 5
lrate = 0.01
decay = lrate/epochs
sgd = SGD(lr=lrate, momentum=0.9, decay=decay, nesterov=False)
model.compile(loss='binary_crossentropy', optimizer=sgd, metrics=['accuracy'])
print(model.summary())
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_1 (Conv2D)            (None, 100, 1104, 1104)   200       
_________________________________________________________________
dropout_1 (Dropout)          (None, 100, 1104, 1104)   0         
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 100, 1104, 1104)   10100     
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 100, 552, 552)     0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 30470400)          0         
_________________________________________________________________
dense_1 (Dense)              (None, 10)                304704010 
_________________________________________________________________
dropout_2 (Dropout)          (None, 10)                0         
_________________________________________________________________
dense_2 (Dense)              (None, 2)                 22        
=================================================================
Total params: 304,714,332
Trainable params: 304,714,332
Non-trainable params: 0
_________________________________________________________________
None

and start training the model:

In [18]:
# Fit the model
model.fit(X_image_train, y_image_train, validation_data=(X_image_test, y_image_test), 
          epochs=epochs, batch_size=1)
# Final evaluation of the model
scores = model.evaluate(X_image_test, y_image_test, verbose=0)
print("Accuracy: %.2f%%" % (scores[1]*100))
Train on 4 samples, validate on 2 samples
Epoch 1/5
4/4 [==============================] - 136s 34s/step - loss: 7.2831 - acc: 0.5000 - val_loss: 8.0590 - val_acc: 0.5000
Epoch 2/5
4/4 [==============================] - 127s 32s/step - loss: 8.0590 - acc: 0.5000 - val_loss: 8.0590 - val_acc: 0.5000
Epoch 3/5
4/4 [==============================] - 127s 32s/step - loss: 4.0295 - acc: 0.7500 - val_loss: 8.0590 - val_acc: 0.5000
Epoch 4/5
4/4 [==============================] - 133s 33s/step - loss: 12.0886 - acc: 0.2500 - val_loss: 8.0590 - val_acc: 0.5000
Epoch 5/5
4/4 [==============================] - 133s 33s/step - loss: 12.0886 - acc: 0.2500 - val_loss: 8.0590 - val_acc: 0.5000
Accuracy: 50.00%

Very good news here is that the CNN undertands the input and does not throw any errors. The output of course does not make sense since we generated our data sets in a strange way, i.e. assigning random classes to random images. This was done on purpose, here now most important is to check that the input can be understood by the model, i.e. there are no syntax errors, so the mission is accomplished.

Case-Control Cell Image Classification with CNN

Here we will implement a simple Case-Control classification with CNN using the augmented data set, where "cases" will be images from the first photograph and "controls" will be augmented images from the second photograph. we will start with reading images, converting them to 4-dimensional numpy arrays and building X_image and y_image data set, where the former contains the numpy arrays of pixels, while the latter contains labels.

In [108]:
import os
import numpy as np
from PIL import Image
tiff_augmented_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_augmented_images/')
X_image = np.zeros(shape = (0,1,100,100))
y_image = list()
for filename in tiff_augmented_files:
    #print('Working with file {}'.format(filename))
    img = Image.open(('/home/nikolay/WABI/S_Aits/tiff_augmented_images/' + filename))
    imgarray = np.array(img)
    imgarray = np.expand_dims(np.expand_dims(imgarray, axis=0), axis=0)
    #print(imgarray.shape)
    X_image = np.concatenate((X_image, imgarray), axis = 0)
    if('f00' in filename):
        y_image.append(0)
    else:
        y_image.append(1)
print(X_image.shape)
y_image = np.array(y_image)
print(y_image.shape)
print(y_image[0:10])
(1080, 1, 100, 100)
(1080,)
[1 1 1 1 0 0 0 0 1 1]

Let us now split the data set into train and test sub-sets using the scikit-learn handy function:

In [110]:
from sklearn.model_selection import train_test_split
X_image_train,X_image_test, y_image_train,y_image_test = train_test_split(X_image, y_image, test_size = 0.2)
print("The dimensions of X_image_train and y_image_train are",X_image_train.shape,"\b and",
      y_image_train.shape,"\b, respectively")
print("The dimensions of X_image_test and y_image_test are",
      X_image_test.shape,"\b and",y_image_test.shape,"\b, respectively")
The dimensions of X_image_train and y_image_train are (864, 1, 100, 100)  and (864,) , respectively
The dimensions of X_image_test and y_image_test are (216, 1, 100, 100)  and (216,) , respectively

Now let us build a simple CNN model where we will use data normalization (by max value) and one-hot encoding of the labels:

In [111]:
# Simple CNN model for CIFAR-10
import numpy
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.layers import Flatten
from keras.constraints import maxnorm
from keras.optimizers import SGD
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.utils import np_utils
from keras import backend as K
K.set_image_dim_ordering('th')

# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)

# normalize inputs from 0-255 to 0.0-1.0
X_image_train = X_image_train.astype('float32')
X_image_test = X_image_test.astype('float32')
X_image_train = X_image_train / 255.0
X_image_test = X_image_test / 255.0

# one hot encode outputs
y_image_train = np_utils.to_categorical(y_image_train)
y_image_test = np_utils.to_categorical(y_image_test)
num_classes = y_image_test.shape[1]
In [119]:
# Create the model
model = Sequential()
model.add(Conv2D(100, kernel_size=(3, 3), input_shape=(1, 100, 100), padding='same', activation='relu', 
                 kernel_constraint=maxnorm(3)))
model.add(Dropout(0.2))
model.add(Conv2D(100, kernel_size=(3, 3), activation='relu', padding='same', kernel_constraint=maxnorm(3)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(512, activation = 'relu', kernel_constraint = maxnorm(3)))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation = 'sigmoid'))
# Compile model
epochs = 5
lrate = 0.01
decay = lrate / epochs
sgd = SGD(lr = lrate, momentum = 0.9, decay = decay, nesterov = False)
model.compile(loss = 'binary_crossentropy', optimizer = sgd, metrics = ['accuracy'])
print(model.summary())
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_3 (Conv2D)            (None, 100, 100, 100)     1000      
_________________________________________________________________
dropout_3 (Dropout)          (None, 100, 100, 100)     0         
_________________________________________________________________
conv2d_4 (Conv2D)            (None, 100, 100, 100)     90100     
_________________________________________________________________
max_pooling2d_2 (MaxPooling2 (None, 100, 50, 50)       0         
_________________________________________________________________
flatten_2 (Flatten)          (None, 250000)            0         
_________________________________________________________________
dense_3 (Dense)              (None, 512)               128000512 
_________________________________________________________________
dropout_4 (Dropout)          (None, 512)               0         
_________________________________________________________________
dense_4 (Dense)              (None, 2)                 1026      
=================================================================
Total params: 128,092,638
Trainable params: 128,092,638
Non-trainable params: 0
_________________________________________________________________
None
In [120]:
import pydot_ng as pydot
import graphviz
from IPython.display import SVG
from keras.utils.vis_utils import model_to_dot
SVG(model_to_dot(model).create(prog='dot', format='svg'))
Out[120]:
G 140544743451560 conv2d_3_input: InputLayer 140544743450888 conv2d_3: Conv2D 140544743451560->140544743450888 140544743036352 dropout_3: Dropout 140544743450888->140544743036352 140544743452008 conv2d_4: Conv2D 140544743036352->140544743452008 140544743037528 max_pooling2d_2: MaxPooling2D 140544743452008->140544743037528 140544743037584 flatten_2: Flatten 140544743037528->140544743037584 140572946602696 dense_3: Dense 140544743037584->140572946602696 140544743456720 dropout_4: Dropout 140572946602696->140544743456720 140544743492576 dense_4: Dense 140544743456720->140544743492576
In [122]:
# Fit the model
estimator = model.fit(X_image_train, y_image_train, validation_split = 0.2, shuffle = True, 
                      epochs = epochs, batch_size = 5)
print("Training Accuracy: ",estimator.history['acc'][-1])
print("Validation Accuracy: ",estimator.history['val_acc'][-1])

# summarize history for accuracy
plt.plot(estimator.history['acc'])
plt.plot(estimator.history['val_acc'])
plt.title('Model Accuracy', fontsize = 20)
plt.ylabel('Accuracy', fontsize = 20)
plt.xlabel('Epoch', fontsize = 20)
plt.legend(['Train', 'Validate'], loc='upper left', fontsize = 20)
plt.show()
# summarize history for loss
plt.plot(estimator.history['loss'])
plt.plot(estimator.history['val_loss'])
plt.title('Model Loss', fontsize = 20)
plt.ylabel('Loss', fontsize = 20)
plt.xlabel('Epoch', fontsize = 20)
plt.legend(['Train', 'Validate'], loc='upper left', fontsize = 20)
plt.show()
Train on 691 samples, validate on 173 samples
Epoch 1/5
691/691 [==============================] - 553s 800ms/step - loss: 0.4644 - acc: 0.7366 - val_loss: 0.2286 - val_acc: 0.8671
Epoch 2/5
691/691 [==============================] - 547s 792ms/step - loss: 0.1015 - acc: 0.9689 - val_loss: 0.0256 - val_acc: 1.0000
Epoch 3/5
691/691 [==============================] - 580s 840ms/step - loss: 0.0377 - acc: 0.9899 - val_loss: 0.0130 - val_acc: 1.0000
Epoch 4/5
691/691 [==============================] - 545s 788ms/step - loss: 0.0150 - acc: 0.9971 - val_loss: 0.0504 - val_acc: 0.9913
Epoch 5/5
691/691 [==============================] - 541s 782ms/step - loss: 0.0150 - acc: 0.9986 - val_loss: 0.0122 - val_acc: 0.9971
Training Accuracy:  0.998552821652
Validation Accuracy:  0.997109825901

We observe a very fast (after only 5 epochs) decrease of the loss function and the accuracy quickly reaches 100% on the validation data set.

In [123]:
# Final evaluation of the model
scores = model.evaluate(X_image_test, y_image_test, verbose=0)
print("Accuracy: %.2f%%" % (scores[1]*100))
Accuracy: 99.54%
In [125]:
ann_pred_probs = model.predict(X_image_test)
In [127]:
ann_pred_probs[0:5]
Out[127]:
array([[  5.40934798e-06,   9.99993205e-01],
       [  1.09973818e-18,   1.00000000e+00],
       [  1.02528960e-04,   9.99935508e-01],
       [  1.00000000e+00,   3.16175092e-13],
       [  8.35562275e-07,   9.99999046e-01]], dtype=float32)
In [129]:
y_image_test[0:5]
Out[129]:
array([[ 0.,  1.],
       [ 0.,  1.],
       [ 0.,  1.],
       [ 1.,  0.],
       [ 0.,  1.]])

Thus we reach a very high accuracy of discriminating between first and second photographs, which is to be expected since those two photograps visually look very different, still good to see the demonstration of Computer Vision. We can see that the predicted labels of the images almost always agree with the true labels.

Object Detection

Images received from Sonja Aits were in C01 format which is an intrinsic Cell Profiler (Cellomics) format. We did image processing and augmentation using Pillow Python module. Now we will implement OpenCV Python module in order to detect cells on the images.

In [2]:
import cv2 
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
plt.imshow(img)
plt.show()

We are going to do some simple image manipulation: turn the image to grayscale, binarize and dilate using custom kernels.

In [3]:
fig = plt.figure(figsize=(20, 15))

#grayscale 
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) 
plt.subplot(131)
plt.imshow(gray)
plt.title('Grayscale')

#binary 
ret,thresh = cv2.threshold(gray,127,255,0) 
plt.subplot(132)
plt.imshow(thresh)
plt.title('Binary')

#dilation 
kernel = np.ones((1,1), np.uint8) 
img_dilation = cv2.dilate(thresh, kernel, iterations=1) 
plt.subplot(133)
plt.imshow(img_dilation)
plt.title('Dilated')
plt.show()

Next, we will find contours using the OpenCV function findCoutours:

In [4]:
# find contours
ctrs, hier = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
# sort contours
sorted_ctrs = sorted(ctrs, key=lambda ctr: cv2.boundingRect(ctr)[0])
In [5]:
cv2.drawContours(img, ctrs, -1, (0,255,0), 3)
Out[5]:
array([[[ 7,  7,  7],
        [ 9,  9,  9],
        [ 6,  6,  6],
        ..., 
        [ 4,  4,  4],
        [ 4,  4,  4],
        [ 6,  6,  6]],

       [[ 9,  9,  9],
        [ 8,  8,  8],
        [ 7,  7,  7],
        ..., 
        [ 6,  6,  6],
        [ 6,  6,  6],
        [ 5,  5,  5]],

       [[ 8,  8,  8],
        [10, 10, 10],
        [ 7,  7,  7],
        ..., 
        [ 4,  4,  4],
        [ 6,  6,  6],
        [ 8,  8,  8]],

       ..., 
       [[32, 32, 32],
        [31, 31, 31],
        [29, 29, 29],
        ..., 
        [35, 35, 35],
        [34, 34, 34],
        [39, 39, 39]],

       [[29, 29, 29],
        [31, 31, 31],
        [26, 26, 26],
        ..., 
        [33, 33, 33],
        [36, 36, 36],
        [43, 43, 43]],

       [[22, 22, 22],
        [30, 30, 30],
        [26, 26, 26],
        ..., 
        [30, 30, 30],
        [37, 37, 37],
        [39, 39, 39]]], dtype=uint8)
In [6]:
for i, ctr in enumerate(sorted_ctrs):
    # Get bounding box
    x, y, w, h = cv2.boundingRect(ctr)

    # Getting ROI
    roi = img[y:y + h, x:x + w]

    # show ROI
    # cv2.imshow('segment no:'+str(i),roi)
    cv2.rectangle(img, (x, y), (x + w, y + h), (0, 255, 0), 2)

    if w > 15 and h > 15:
        cv2.imwrite('/home/nikolay/WABI/S_Aits/bounding_boxes/i.png'.format(i), roi)
fig = plt.figure(figsize=(20, 15))
plt.imshow(img)
plt.show()

Now we have detected the cell on the images, put bounding boxes around them, croped the cell amd written them to separate imge files. Next step is to annotate the bounding boxes and input them to Faster RCNN.

In [ ]: